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ABSTRACT 

A procedure that models coupled thermo -mechanical deformations of viscoelastic rubber 
cylinders by employing the ABAQUS finite element code is described. Computational 
simulations of hysteretic heating are presented for several tall and short rubber cylinders both 
with and without a steel disk at their centers. The cylinders are compressed axially and are then 
cyclically loaded about the compressed state. The non-uniform hysteretic heating of the rubber 
cylinders containing a steel disk is presented. The analyses performed suggest that the coupling 
procedure should be considered for further development as a design tool for rubber degradation 
studies. 
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1. INTRODUCTION 

1.1 Background Rubber is employed to carry large loads in tires, gaskets, and tank track pads. 
It is also used to provide damping and system stability in complex mechanical systems such as 
helicopter rotors. In these applications, the rubber is typically st if fened by the addition of carbon 
black. The filled rubber tends to be a poor conductor of heat, yet it also exhibits very large 
hysteretic energy loss during cyclic loading. Also, the mechanical properties of rubber are 
strongly dependent on temperature. Faced with the above issues, designers interested in 
modeling the detailed response of complex- shaped rubber components need to be able to 
compute the coupled thermo -mechanical behavior of rubber. 

An example of the importance of understanding the thermo -mechanical response of fil led rubber 
is given in a series of papers presented at the "Thirty Second Sagamore Army Materials Research 
Conference" held a Lake Luzerne, NY in 1985 (1-4). In these papers, hysteretic heating, thermo- 
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mechanical degradation, and fatigue of rubber- coated road wheels and tank track pads are all 
discussed. Uncoupled thermo -mechanical finite element analyses, and sensitivity studies were 
conducted with finite element codes. It was observed that the viscoelastic properties and the 
shape of the rubber solid are the most important factors in determining temperature rise (1). The 
degradation studies indicate that the failure of cyclically loaded "rubber- like" polyurethane 
blocks depends on the hard segment transition temperature (2). Experiments were conducted 
which proved that the large strain hysteretic heating rate does not correlate with the heating rates 
predicted using the popular complex modulus material data (3). Also determined is the fact that 
failure under cyclic loading can be "significantly different from that obtained in constant rate 
testing" (4). These conclusions suggest that detailed computational simulations of large strain 
dynamic loading of rubber-like solids, performed as part of a material degradation study, require 
accurate modeling of the large strain viscoelastic properties and a coupling of the mechanical and 
thermal models. 

The purpose of this work is to establish and test a computational tool for analyzing hysteretic 
heating in rubber components. In this paper, large deformation rubber viscoelasticity and heat 
transfer finite element models are coupled and employed to computationally simulate the 
hysteretic heating of dynamically loaded rubber cylinders. 

The formulations for large strain rubber viscoelasticity and heat transfer employed in ABAQUS 
are outlined below to facilitate the description of the thermo -mechanical coupling performed in 
this study. Detailed information on these two formulations and their finite element 

implementation is available in the ABAQUS Theory Manual (5). Additional information on 

formulations for rubber viscoelasticity are available in the literature (6-10). Only moderate 
temperature changes were simulated in this study, so time -temperature superposition is not 
discussed. 

1.2 Rubber Viscoelasticity The virtual work statement, without inertial effects included, for a 
solid of volume V and surface area S is: 

8Wj = js :SD V dV = jSv r t dS + JSv r f dV [1] 

V S V 

where 5 W, is the internal energy due to the virtual displacement Sv , s is the Cauchy stress, 
8D V is the virtual rate of deformation, t is the traction stress vector acting on S , and f is the 

body force vector acting within V . Details on the notation used below are given in the 
Appendix. Equation [1] is used to build the mechanical finite element equilibrium equations. 
The following details describe how the viscoelastic behavior of rubber is approximated. In the 
reference configuration, the strain energy part of Equation [1] is written: 

bW, = J/s :SD V dV () [2] 

v 0 

where V 0 is the volume V in its reference state. When the solid is rubber. Equation [2] is 
expressed as: 
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where U is the solid's strain energy density function. In this effort, a standard polynomial form 
of the rubber's strain energy density was employed as follows: 
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where the constants C, ; and D ; are determined from measured stress- strain data, and e th is the 

thermal expansion strain. A history integral method is employed in ABAQUS to model 
viscoelastic material behavior. In the viscoelasticity formulation the material constants in the 

strain energy density function are made time dependent as follows: 
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where the constants g e , k,,, and x e are determined by fitting the analytical form of relaxation 

stresses to experimental data (5, 9-12), the constants C \ ] and D ( - are the "elastic" constants, 

mention above in Equation [4], determined from high strain rate loading tests, and the superscript 
0 refers to the material's instantaneous response. 


1.3 Heat Transfer. The variational statement of the energy balance equation for heat transfer, 
together with Fourier's law, for a deformed solid of volume V , and surface area S is: 


{P 


^480 dV+ f^-k— dV = [80 rdV+\8QqdS 
dt J dx dx J ] 


[ 6 ] 


where 0 is the temperature, U e is the internal thermal energy, p is the density, q is the heat 

flux per unit area, r is the heat supplied per unit volume, k is a conductivity matrix, and 80 is a 
virtual temperature field satisfying the essential boundary conditions. Equation [6] is used to 
build the transient heat transfer finite element equations. 


2. COUPLED THERMO-MECHANICAL MODEL 


Finite element discretization of the mechanical and thermal variational statements, Equations [1] 
and [6], results in systems of time dependent mechanical and thermal differential equations. For 
the variational statements given above, the only thermo -mechanical coupling that exists is 
through the thermal expansion strain, £ f/i . Additional thermo -mechanical coupling can be 

approximated by computing the rate of viscoelastic energy dissipation from the mechanical 
equations and inputting that rate into the thermal equations. ABAQUS internally computes and 
stores a running total of the energy dissipated as a function of time. Internal variables can be 
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employed to compute the energy dissipation rate from this running total and pass it on to the 
thermal analysis. 

2.1 Element Hie eight- node hybrid axisymmetric CAX8RHT element in ABAQUS was uti liz ed 
for the examples given below. The element uses biquadratic displacement interpolation, and 
bilinear temperature and pressure interpolation. It also employs reduced integration (see Figure 
1 .) 



r 

Figure 1. Axisymmetric CAX8RHT element. 


The variables employed in the thermal element are the temperatures, 9 , at the four comer nodes. 
The stress element has radial, u r , and axial, u z , displacements at the eight nodes and hydrostatic 
pressure, p l , variables at each of the four integration points. Heating rates are applied at the four 
integration points. 

2.2 Coupling Procedure Two internal variables are introduced at each of the four integration 
points in the element. When convergence is achieved, at the end of a time increment, the energy 
dissipated during the time increment is computed internally (for each integration point in each 
element) by ABAQUS. The increment is added to and stored as a running total at each 
integration point. The "user subroutine" option in ABAQUS was employed with the two internal 
variables to compute the energy dissipation rate and input the rate to the thermal analysis. At 
each integration point, one internal variable stores the "running total" dissipated energy at the 
beginning of the time interval. These stored values are used with the ABAQUS "running totals" 
at the end of the time interval to compute the rates of energy dissipation across the time interval 
at each integration point. The rate of energy dissipation is stored as the second internal variable 
and is sent forward to the next time interval to serve as a heating rate in the thermal analysis. 

2.3 Cylinder Dimensions, Material Properties, and Loading To investigate the application of 
the procedure described above, eight cylinders were analyzed for hysteretic heating. Only 
moderate temperature changes were being simulated so the elastic material constants were not 
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treated as temperature dependent. Having the elastic constants independent of temperature 
simplified the calculations presented in this effort, but it is not a lim itation of the algorithm. 

2.3.1 Dimensions There were two groups of four cylinders each. One group consisted of 
uniform cylinders and the other group had steel disks at their centers, see Figure 2. All cylinders 
had the radius, R = 0.0282 m. There were four cylinder heights in each group. The heights were 
H = 0.05 m, 0.0375 m, 0.025 m, 0.0125 m respectively. The cylinders were compressed between 

steel fixtures. The model simulates the case when a lubricant maintains the fixture- rubber 
interface as ftictionless. The internal steel disks were completely attached (bonded) to the 
rubber. The height and radius of the disks were 0.0025 m , and 0.0141 m , respectively. 


Axisymmetric FE Model 

Internal heat generated 
by the dissipated 
viscoelastic energy 

Steel disk 

Rubber cylinder _____ 

Heat convects to air 
at outer surface 

Heat conducts to fixture - 
at interface 

Figure 2. Rubber cylinder, finite element mesh, fixture, and steel disk. 

2.3.2 Rubber Properties The rubber energy density was modeled with a two-term Mooney- 
Rivlin (6) expansion and with two terms in each Prony series (see Equations [4] and [5]). The 
material constants employed are representative of a soft rubber and are listed below (13)- (15). 

The viscoelastic constants are: 

Cl o = 454545.45 Pa, = 45454.54 Pa, and D,° = 6.667 *10“ 10 Pa~ l 
g, = 0.2, k { = 0.2, x j = 0.2 s 
and g 2 =0.2, k 2 =0.2,t 2 =1.05 

The thermal material constants are: 

Conductivity K =0.20934 J/[°Cm s) 

Density p = 1000. kg/ m 3 

Specific heat c = 2093 .4 J /[kg °c) 

Expansion a - 80. * 10 -6 (° c) 
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The fil m heat transfer coefficients for the rubber- air and rubber- steel interfaces are 
h A =5.442847 / (c m s) and h s =20934 J /{°C ms), respectively. 

2.3.3 Steel Properties 

The elastic constants are: 

Young’s modulus E = 206.8 GPa 

Poisson’s Ratio v = 0.3 


The thermal constants are: 
Conductivity 
Density 
Specific heat 
Expansion 


K =45.83379 J/(°Cm sec) 
p =7849. kg/m 3 
c = 460. J/(kg°c) 
a = 12. * 10 _6 (°c) 


2.3.4 Loading Each cylinder was subjected to a constant compressive load with a superimposed 
cyclic load of sufficient magnitude to produce large strain hysteresis. The loading was: 

f{t) = -2100 -1400 sin (40.84 1) N [7] 

where t is the time in seconds. The maximum deformation of the tall ( H - 0.05 m ) uniform 
cylinder after 130 cycles is shown in Figure 3. The mesh refinement near the top and outer 
edges is to accommodate the heat transfer gradients at those locations. 



The loading and the displacement of the top of the cylinder are shown in Figure 4 for a time 
interval of 1.0 s. As expected, the displacement curve demonstrates viscoelastic softening or 
"cyclic creep." 
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Figure 4. Viscoelastic softening of a tall uniform cylinder (H = 0.05 m). 

3. COMPUTED TEMPERATURE DISTRIBUTIONS 

A number of analyses were performed to evaluate the nonlinear elastic, the viscoelastic, and the 
transient heat transfer finite element algorithms separately. Hand calculations and finite 
difference calculations verified that the algorithms functioned correctly. In the prehminary 
calculations, the coupling procedure was used with properties representative of track pad 
materials and the rate at which the temperature increased was similar to rates measured for tests 
on track pad materials (2). The thermal boundary conditions did not significantly affect the 
temperature fields computed. The following analyses were performed, employing the coupling 
procedure described above, to investigate the heating of soft rubber cylinders undergoing large 
strain dynamic deformations. 

3.1 Uniform Cylinders The cyclic load given by Equation [7] was applied to the four uniform 
cylinders described above. The temperatures, at the center of the cylinders as a function of time 
for the first 20 s of loading are presented in Figure 5. The frictionless mbber- fixture surface 
allows the strains to be uniform in the cylinder. Rubber is a poor conductor of heat and the fact 
that the cylinders heated nearly uniformly regardless of shape (tall or short) was expected. 

3.2 Cylinders Containing a Disk It is difficult to estimate hysteretic heating in mbber solids of 
complex shape because coupled thermo -mechanical analyses are needed in regions of high strain 
gradients. The group of cylinders containing disks were cyclically loaded to observe the ability 
of the coupling procedure to predict the distribution of viscoelastic heating in a mbber solid with 
high strain gradients. The results obtained are reasonable. Figure 6 shows the meshes on the 
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reference and deformed shapes for the tall cylinder. The temperature distribution in the tall 
block after 20 s of dynamic loading is shown in Figure 7. Temperatures at points A, B, C, and D 
in Figure 7 are plotted as a function of time in Figure 8. The outer radial end of the internal disk 
(point C) is predicted to heat much faster than the rest of the cylinder. Snail oscillations in the 
computed heating rates were observed in these computations. The cause of these oscillations has 
not been determined. 



Figure 5. Temperature as a function of time at the center of each uniform cylinder. 


Steel disk 


Surface film for 
heat conduction 



Frictionless boundary 


Surface film for 

heat convection 


Figure 6. Deformation of a tall cylinder with an internal disk (H = 0.05 m). 
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Figure 7. Temperature distribution in a tall cylinder with an internal disk (H = 0.05 m). 
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Figure 8. Temperature as a function of time for a tall cylinder with an internal disk (H = 0.05m). 



The results were similar for the shorter cylinders. The temperature distribution in the shortest 
cylinder is presented in Figure 9. As in the case of the tall cylinder, the region of high strain 
gradient located near the outer radial end of the internal disk (point C) has the most rapid rise in 
temperature (see Figure 10.) 
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Steel disk 


Figure 9. Temperature distribution in a short cylinder with an internal disk (H = 0.0125 m). 



Time (s) 

Figure 10. Temperature as a function of time for a short cylinder 
with an internal disk (H = 0.0125 m). 


4. CONCLUSIONS 

Accurate predictions of the strain and temperature distributions in rubber components, employed 
in dynamically loaded structures, are required to perform degradation studies. A procedure 
which couples a viscoelastic stress analysis with a heat transfer analysis was described. The 
ABAQUS finite element code was employed to demonstrate the procedure. Both the stress and 
the thermal analyses are valid for large strains and displacements. A user subroutine was written 
to track dissipation of energy across time intervals and determine heating rates. The heating rate 
data is passed forward one time interval in the procedure. 
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Hie thermo -mechanical heating of tall and short uniform rubber cylinders (without internal 
disks) was computed with the coupling procedure. The viscoelastic material properties 
employed are valid for large strain deformations of soft rubber. The time dependent strains in 
the cylinders were uniform, and uniform heating was computed. 

Hie thermo -mechanical heating of tall and short rubber cylinders with internal steel disks was 
also studied. Hie internal steel disks provided high strain gradients within the rubber cylinders 
and nonuniform hysteretic heating is observed. 

The analyses performed suggest that the coupling procedure should be considered for further 
development as a design tool for rubber degradation studies. An integrated effort in which the 
large strain viscoelastic material properties, the heat transfer properties, and the heating rates due 
to cyclic loading are all determined is recommended. 
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APPENDIX 


The following notation is included to assist the reader with the description of the ABAQUS finite 
element algorithm for rubber viscoelasticity. 


Xf (i =1,2,3) 
x t (i = 1,2,3) 
x = x(X,t) 

ax 


coordinates of a material point in the reference configuration, 
coordinates of a material point in the deformed configuration, 
vector mapping between the reference and deformed configurations. 

deformation gradient. 


J= det(F) 

__L 

F = / 3 F 
B = FF r 

I x = /t(b) 

1=^((1 2 Mbb)) 


determinate of F which measures volume change. 

deformation gradient scaled for volume change. 

left Cauchy Green strain tensor. 

first strain invariant (adjusted for volume). 

second strain invariant (adjusted for volume). 


8u 


SL = 


dSu 

dx 


displacement, 
gradient of displacement. 


SD = ^(SL + 5L r ) 


rate of deformation. 


SD V 

A :B 

5e vo/ = 1 : 5D 


rate of deformation computed using virtual displacement 5v . 
scalar product of matrices A and B . 
volumetric strain rate. 


Se =SD-— Se vo, I 
3 


P = 



: s 


deviatoric strain rate, 
pressure stress (hydrostatic). 
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